The present invention relates to exploration seismic reflection surveying and more particularly, relates to the processing of exploration seismic reflection data to enhance information in seismic signals reflected from contrasts or differences in elastic constants or densities in the subsurface of the earth.
The methods of the present invention which are described herein are generally discussed in terms of compressional-wave (PP) seismic data acquisition and processing, which is the most common form of seismic data used in exploration seismology. However, it should be understood that these methods are equally applicable to shear-wave seismic data.
Conventional land or marine seismic acquisition techniques involve the use of an appropriate source to generate seismic energy and a set of receivers, spread out along or near the surface of the earth on land or at or near the water surface or water bottom in a water covered area, to detect any reflected seismic signals due to seismic energy stricking subsurface geologic boundaries. These signals are recorded as a function of time and subsequent processing of these time varying signals, i.e. seismic "traces" or seismic data, is designed to reconstruct an appropriate image of the geologic boundaries of the subsurface and to obtain information about the subsurface materials. In simplistic terms, this conventional process has a seismic wave, from a source of seismic energy, travelling down into the earth, reflecting from a particular geologic interface (i.e. a change or contrast in elastic constants and/or densities), and returning to the surface, where it may be detected by an appropriate receiver.
If the seismic-wave velocity is known as a function of depth and lateral position, and if the position and dip of a plane geologic interface are known, the time for the wave to travel down to that particular reflecting interface and reflect back to the surface can be computed for any source and receiver locations. This two-way travel time are usually described by a function t(X,Z), where Z is the depth to the reflecting interface (contrast in elastic constants or density) and X is the horizontal distance (offset) between source and receiver.
If the elastic constants and densities of the materials above and below a planar reflecting interface are known, then the reflection coefficient may be computed. This reflection coefficient is the ratio of reflected amplitude to incident amplitude and will depend on the angle of incidence at the reflecting interface. The angle of incidence, .theta., is the angle between the ray normal to the incident downgoing wavefront and a line normal to the interface; as is well known, the incident and reflected rays will be in a plane normal to the interface. This angle of incidence increases with increasing offset X. The reflection coefficient for a compressional wave from a particular interface will be designated by the function R.sub.P (.theta.), where the angle .theta. may be related to the offset distance X and depth of reflector Z by raytracing if the compressional-wave velocity at all points in the earth is known; this velocity information, or a reasonable approximation thereto, is referred to as a "velocity model". For a given reflector, the reflection angle, .theta., and offset, X, are geometrically related, so any discussions herein in terms of dependence upon offset (offset dependence) is equivalent to dependence upon reflection angle (angular dependence). The angular (or offset) dependence of reflection amplitude may be computed exactly for a point source and plane reflector, however in most practical cases it may be approximated adequately by plane-wave reflection coefficients (reflection coefficients for an incident plane wave) which are easily calculated using expressions derived from the results of Zoeppritz (see for example, Bull Seismol. Soc. Am., Vol. 66, 1976, pp. 1881-1885, Young, G. B. and Braile, L. W.). For a compressional-wave reflection from a planar interface between two media having a small contrast (i.e., with the medium containing the incident and reflected waves having a compressional velocity V.sub.P, a shear velocity V.sub.S, and a density .rho., and the other medium having a compressional velocity of V.sub.P +dV.sub.P, a shear velocity V.sub.S +dV.sub.S, and a density .rho.+d.rho., and where dV.sub.P /V.sub.P, dV.sub.S /V.sub.S, and d.rho./.rho. are small compared to one), the offset (or reflection angle) dependence of reflection amplitude may be described for angles of incidence less than the critical angle by an expansion of the form, EQU R.sub.P (.theta.)=R.sub.P (O)+K sin.sup.2 (.theta.)+L sin.sup.4 (.theta.)+ . . . (1)
For the discussion herein, the angles of incidence are limited to angles such that the terms of the order of sin.sup.4 (.theta.) and higher are negligible. In Equation (1), R.sub.P (O) is the normal incidence (.theta.=O) reflection coefficient; R.sub.P (O) depends only on the densities and compressional velocities of the two media. K is a constant, which also depends on the elastic properties and densities of the media. The relationship of K to the elastic properties and densities may be expressed in a number of ways. One particularly simple expression which relates K to the contrasts in shear velocities, compressional velocities, and densities is EQU K=R.sub..alpha. -4(V.sub.S /V.sub.P).sup.2 (2R.sub..beta. +R.sub..rho.),(2)
where EQU R.sub..alpha. =dV.sub.P /(2V.sub.P +dV.sub.P), (2a) EQU R.sub..beta. =dV.sub.S /(2V.sub.S +dV.sub.S), and (2b) EQU R.sub..rho. =d.rho./(2.rho.+d.rho.). (2c)
Also, in terms of these same coefficients, the normal incidence or zero offset reflection coefficients are given exactly by, ##EQU1## For sufficiently small values of R.sub..alpha., R.sub..beta., and R.sub..rho., equations 2d and 2e may be approximated as, EQU R.sub.P (O)=R.sub..alpha. +R.sub..rho., and (2f) EQU R.sub.S (O)=R.sub..beta. +R.sub..rho.. (2g)
Thus, measurement of the normal incidence compressional-wave reflection coefficient, R.sub.P (O), gives information about the densities and compressional velocities, while measurement of the offset dependence constant K can provide information about the densities and shear velocities of the media. The formulas given above are for small contrasts in the elastic properties and densities above and below the planar interface; of course, more general theoretical relations may be used. Similar relationships are well known for the offset dependence of shear-wave reflection coefficients, although the particular form for such shear-wave equations that are analogous to Equation 2 is quite different.
There are a number of geologic questions important to exploration for hydrocarbons which can be answered by acquiring a knowledge of both the compressional- and shear-wave properties of the subsurface materials. For instance, these materials are generally porous with various fluids filling the pore space. The compressional velocity of a seismic wave in such media depends strongly on the rock matrix properties as well as those of the pore fluid. On the other hand, shear-wave velocities depend strongly on the rock matrix but only slightly on the pore fluid. Thus, detailed study of both the compressional and shear properties of the media provides an opportunity to characterize any changes in seismic response as being due to changes in fluid content (e.g. from brine to oil, or oil to gas) or changes in the rock matrix (e.g. from sandstone to shale or a change in porosity). The ratio of V.sub.P to V.sub.S is often a useful diagnostic feature of such changes. It should be noted that, even without lateral variation, in many cases the recognition of fluid content or rock type may be possible with an accurate knowledge of the compressional and shear properties at a single location. Distinguishing between fluid effects and lithology effects, and detecting different porosity and lithology types are of vital exploration interest and the desire to make such distinctions has engendered significant effort in the measurement and interpretation of shear properties in addition to the information concerning compressional properties traditionally inferred from conventional PP reflection prospecting.
It is generally the objective of seismic exploration to generate seismic energy, make measurements of the reflection amplitude of this energy at various offsets and for various times, and then, by employing various processing steps on this seismic data, to deduce the geometry as well as some of the elastic properties and densities of the materials of the earth through which the seismic energy has propagated and from which it has been reflected.
One such processing method which results in some knowledge of a velocity model (that is the velocity of seismic waves as a function of Z) is conventional moveout velocity analysis. This velocity model may then be used for dynamic correction followed by stacking the corrected seismic data, as described below.
A seismic reflection from an interface will arrive at a receiver after a two-way travel time, denoted and used herein as, t(X), where X is defined and used herein as the distance between the source and the receiver, or "offset" distance. This "moveout time" t(X) may be used to "dynamically correct" seismic data acquired at an offset distance X so that a reflection is adjusted in time to appear as if it had been acquired at zero offset, i.e. X=0. Conventional "stacking" is accomplished by summing or averaging such dynamically corrected data.
For seismic data where the downward and upward propagating wave types are the same, either compressional or shear, the square of this moveout time has a well known expansion about X=0: EQU t.sup.2 (X)=t.sub.0.sup.2 +a.sub.2 X.sup.2 +a.sub.4 X.sup.4 + . . .,(3)
where the coefficients a.sub.2, a.sub.4, a.sub.6, etc. . . . are known functions of the propagation velocity at each depth and t.sub.0 is the X=0 (zero offset) two-way travel time. The propagation velocity at a particular depth is known as the "interval" velocity at that-depth. Experience with compressional-wave data has indicated that in most areas the P-wave vertical velocity distributions are such that, within the offset range (X&lt;Z) and frequency range (5-100 Hz) normally considered in exploration seismology, terms of higher power than X.sup.2 are small and for many purposes may be ignored. Ignoring the higher order terms is equivalent to approximating the stratified earth with a single horizontal layer having an "effective" velocity, V.sub.e. For this simplified case Equation 3 reduces to the "normal hyperbolic moveout" expression EQU t.sup.2 (X)=t.sub.0.sup.2 +X.sup.2 /V.sub.e.sup.2, (4)
where ##EQU2## and EQU t.sub.0 =2.intg.[dZ/V(Z)]. (4b)
For the simplistic discussions hereinabove, the earth model is assumed to be one for which velocities may vary with depth, but not laterally; however, the methods of the present invention may be applied in more general cases. For example, this model and analysis is also appropriate for reflections from dipping reflectors as long as the data is gathered by common midpoint (CMP) between source and receiver pairs (as more fully described below), and the degree of dip is not severe, and the reflectors are laterally continuous. If these assumptions are not met, other conventional processing steps, such as migration, may be required, although they are not described herein. For a dipping reflector, the interpretation of the coefficient a.sub.2 of equation 3 as the square of the reciprocal of effective velocity must be modified to include the cosine of the dip angle. For this case equation 4 becomes: ##EQU3## where .alpha. is the dip angle.
Conventional processing of compressional-wave data uses data collected with many sources and many receivers and gathers the data by the common midpoint (CMP) technique, as illustrated in FIG. 1A. Traces with a common "midpoint" between the source and receiver are collected or gathered at a surface gather point. For example, in FIG. 1A, S.sub.1 and R.sub.1 are the source and receiver pair for the first trace and have a midpoint at the surface point (0). FIG. 1B depicts the corresponding hyperbolic moveout of such data (where the numbers used correspond to the subscripts used in FIG. 1A) and FIG. 1C depicts the corresponding variation of reflection coefficient with offset for such a case. That is, after the data have been acquired for a number of sequential source positions the traces for various source-receiver combinations are sorted or gathered into different midpoint groups which have the same or "common" surface location of the "midpoint" between the source and receiver positions (as depicted in FIG. 1A). This sorted or gathered midpoint data is then analyzed or processed (via equations 4) to determine effective velocities V.sub.e, also called normal moveout velocities, for reflections from various depths, i.e. various values of t.sub.o. One method is to determine for each t.sub.o a value of V.sub.e which provides dynamic corrections which maximize the resultant amplitudes of the "stacked" data in a time gate around t.sub.o. Many methods for moveout velocity analysis and dynamic correction of compressional (PP) and shear (SS) seismic data have been in use since the late 1960's, as described for instance by Tanner and Koehler ("Velocity Spectra-Digital Computer Derivation And Application Of Velocity Functions", Geophysics 34, pp. 859-881 (1969) Tanner, M. T. and Koehler, F.).
The original basis for CMP processing is the fact that each trace in a gather images (or consists of reflections from) approximately the same subsurface points, and, when properly adjusted for differing path lengths, the set of corrected traces may be averaged to give an enhanced picture of the reflection response of the earth below that CMP surface location by emphasizing true primary reflections and discriminating against multiple reflections and other undesirable noise. It is usually assumed that the "stacked" trace represents the normal incidence (zero-offset) response of the earth. While this procedure has been very effective in improving signal-to-noise ratios for seismic data in many areas, it ignores the fact that reflection amplitudes vary as a function of offset and that the stacked trace is not equivalent to a normal incidence trace.
More recently, methods have been described for measuring and interpreting the variation with offset of the reflection amplitude from a given subsurface interface. Techniques which are taught for example in U.S. Pat. No. 4,562,558 to Ostrander, 4,573,148 to Herkenhoff et al, 4,570,246 to Herkenhoff et al, 4,316,267 to Ostrander, 4,316,268 to Ostrander, and 4,534,019 to Wiggins et al explicitly describe methods for measuring and interpreting amplitude variation with offset. Since reflections appear at increasing times at increasing offsets, all these newer methods also require a technique for identifying the time at which a given reflection appears at each offset. Various conventional methods have been used for identifying these times in the above-noted patents, as well as in other publications on the subject. One method simply proposes following a particular trace feature, say a maximum amplitude as it moves out with offset. A more common method utilizes formulas like equation (4) with an effective velocity determined by conventional means (i.e. ignoring amplitude versus offset) as described in the reference noted hereinabove.
Applicants have discovered that in many cases neither the conventional technique of event moveout determination by conventional velocity analysis, or the conventional technique of following the peak amplitude (or other feature) of an event as a function of offset is of sufficient accuracy for use in studying reflection coefficient variation with offset, specifically because these techniques do not properly account for the effects on moveout determination of amplitude variation with offset. Examples, where these conventional techniques fail, include situations where the reflection amplitude varies with offset in such an extreme manner that the reflection amplitude changes sign, going from positive values to zero and then to negative values over the ranges of offsets being considered. Another case is described later herein, where the varying offset dependence of closely spaced reflectors can have the appearance of a single reflector convolved with a wavelet which changes with offset in a spurious manner. Thus, a method for determining the effective velocity of seismic events more accurately by specifically including the effects of amplitude variations with offset is needed. Moreover, such a method could simultaneously determine a more accurate measure of amplitude variation with offset and normal incidence reflection amplitude.
These and other limitations and disadvantages of the prior art are overcome by the present invention, however, and improved methods are provided for determining from seismic data normal moveout velocities and associated dynamic corrections that preserve amplitude versus offset information to provide better estimates of effective velocities, the normal incidence reflection amplitudes, and the offset dependence of the reflection amplitudes.